rm(list=ls())  
if(Sys.info()['user']=='dannychoi'){pathData<-"/Users/dannychoi/Dropbox/Projects/DCMPNS_Language/DataAnalysis"}
if(Sys.info()['user']=='doc43'){pathData<-"/Users/doc43/Dropbox/Projects/DCMPNS_Language/DataAnalysis"}
setwd(pathData)
library(stargazer)
library(estimatr)
library(readr)
library(lubridate)
library(tidyverse)
library(Rmisc)
library(lfe)

# Importing merged data
load(file="merged.RData")

##########################################
##### Main Paper: Figures and Tables #####
##########################################

##### Main Figure 3: Dicrimination Plot

mergeddisttest <- Rmisc::summarySE(merged, measurevar="help", groupvars=c("language"))
mergeddisttest$point <- sprintf(mergeddisttest$help*100, fmt = '%#.2f')

# No hijab vs native comparison
summary(lm(help~nhvn, data=merged))

# Hijab vs native comparison
summary(lm(help~hvn, data=merged))

# Hijab vs no hijab comparison
summary(lm(help~hvnh, data=merged))



# Figure 3
interval1 <- -qnorm((1-0.90)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier
mylabels <- c("(1)\nImmigrant Hijab\n(N=618)", "(2)\nImmigrant No Hijab\n(N=638)", "(3)\nNative\n(N=312)")

library(forcats)
pd <- position_dodge(0.1)
pdf("disttest_iteration_merged.pdf", width = 7, height = 5)
mergeddist <- ggplot(mergeddisttest, aes(x=factor(language), y=100*help, fill=factor(language))) +
  geom_bar(position=pd, stat="identity", width=0.45, alpha=0.9) +
  geom_errorbar(aes(ymin = 100*(help - se*interval2), ymax = 100*(help + se*interval2)),
                width = 0.1, position = position_dodge(width = 1/2)) +
  geom_text(aes(label = point), size = 4.0, vjust=-0.5, hjust = -0.4, nudge_x = 0.05) +
  ggtitle(" ") +
  xlab("") +
  ylab("Assistance Rates (%)") +
  ylim(0.0, 110) +
  scale_x_discrete(breaks=c("1", "2", "3"), labels=mylabels)+
  theme_minimal() +
  theme(axis.text.y=element_text(colour = "gray30", size=11), axis.ticks.y=element_blank()) +
  theme(axis.text.x=element_text(colour = "gray30", size=11), axis.ticks.x=element_line(colour="gray30"))
mergeddist<- mergeddist + geom_segment(aes(x=2, y=87, xend=2.2, yend=90), colour="darkgray") +
  geom_segment(aes(x=2.2, y=90, xend=2.8, yend=90), colour="darkgray") +
  geom_segment(aes(x=2.8, y=90, xend=3, yend=87), colour="darkgray")
mergeddist <- mergeddist + annotate("text", x=2.5, y=93, label="-2.7%p (p=0.347)", size=3)
mergeddist<- mergeddist + geom_segment(aes(x=1, y=83, xend=1.2, yend=86), colour="darkgray") +
  geom_segment(aes(x=1.2, y=86, xend=1.8, yend=86), colour="darkgray") +
  geom_segment(aes(x=1.8, y=86, xend=2, yend=83), colour="darkgray")
mergeddist <- mergeddist + annotate("text", x=1.5, y=89, label="12.2%p (p<0.001)", size=3)
mergeddist<- mergeddist + geom_segment(aes(x=1, y=98, xend=1.2, yend=101), colour="darkgray") +
  geom_segment(aes(x=1.2, y=101, xend=2.8, yend=101), colour="darkgray") +
  geom_segment(aes(x=2.8, y=101, xend=3, yend=98), colour="darkgray")
mergeddist <- mergeddist + annotate("text", x=2, y=104, label="9.5%p (p=0.003)", size=3)
mergeddist <- mergeddist + scale_fill_manual(values=c("#CC0000","#CC33CC", "#0000FF"))
mergeddist <- mergeddist + theme(legend.position="none")
print(mergeddist)
dev.off()



##### Main Figure 4: Language Effects
mergedofftest <- Rmisc::summarySE(merged, measurevar="help", groupvars=c("languageoff"))
mergedofftest$point <- sprintf(mergedofftest$help*100, fmt = '%#.2f')

# Hijab foreign language vs Hijab German language
summary(lm(help~hfvhg, data=merged))

# No hijab foreign language vs No hijab German language
summary(lm(help~nfvng, data=merged))

# Hijab German language vs Native 
summary(lm(help~hgvn, data=merged))

# Hijab German language vs No hijab German language
summary(lm(help~hgvnhg, data=merged))

# Figure 4
interval1 <- -qnorm((1-0.90)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier
mylabels <- c("(1)\nHijab\nForeign Language\n(N=284)", "(2)\nHijab\nGerman Language\n(N=334)", "(3)\nNo Hijab\nForeign Language\n(N=309)",
              "(4)\nNo Hijab\nGerman Language\n(N=329)", "(5)\nNative\nGerman Language\n(N=312)")


pd <- position_dodge(0.1)
pdf("langoffset_iteration_merged.pdf", width = 8, height = 5)
mergedoffset <- ggplot(mergedofftest, aes(x=factor(languageoff), y=100*help, fill=factor(languageoff))) +
  geom_bar(position=pd, stat="identity", width=0.45, alpha=0.9) +
  geom_errorbar(aes(ymin = 100*(help - se*interval2), ymax = 100*(help + se*interval2)),
                width = 0.1, position = position_dodge(width = 1/2)) +
  geom_text(aes(label = point), size = 4.0, vjust=-0.5, hjust = -0.4, nudge_x = 0.05) +
  ggtitle(" ") +
  xlab("") +
  ylab("Assistance Rates (%)") +
  ylim(0.0, 110) +
  scale_x_discrete(breaks=c("1", "2", "3", "4", "5"), labels=mylabels)+
  theme_minimal() +
  theme(axis.text.y=element_text(colour = "gray30", size=11), axis.ticks.y=element_blank()) +
  theme(axis.text.x=element_text(colour = "gray30", size=11), axis.ticks.x=element_line(colour="gray30"))
mergedoffset<- mergedoffset + geom_segment(aes(x=1, y=77, xend=1.2, yend=80), colour="darkgray") +
  geom_segment(aes(x=1.2, y=80, xend=1.8, yend=80), colour="darkgray") +
  geom_segment(aes(x=1.8, y=80, xend=2, yend=77), colour="darkgray")
mergedoffset <- mergedoffset + annotate("text", x=1.5, y=83, label="0.07%p (p=0.984)", size=3)
mergedoffset<- mergedoffset + geom_segment(aes(x=3, y=87, xend=3.2, yend=90), colour="darkgray") +
  geom_segment(aes(x=3.2, y=90, xend=3.8, yend=90), colour="darkgray") +
  geom_segment(aes(x=3.8, y=90, xend=4, yend=87), colour="darkgray")
mergedoffset <- mergedoffset + annotate("text", x=3.5, y=93, label="3.3%p (p=0.320)", size=3)
mergedoffset<- mergedoffset + geom_segment(aes(x=2, y=95, xend=2.2, yend=98), colour="darkgray") +
  geom_segment(aes(x=2.2, y=98, xend=3.8, yend=98), colour="darkgray") +
  geom_segment(aes(x=3.8, y=98, xend=4, yend=95), colour="darkgray")
mergedoffset <- mergedoffset + annotate("text", x=3, y=101, label="13.7%p (p<0.0001)", size=3)
mergedoffset<- mergedoffset + geom_segment(aes(x=2, y=103, xend=2.2, yend=106), colour="darkgray") +
  geom_segment(aes(x=2.2, y=106, xend=4.8, yend=106), colour="darkgray") +
  geom_segment(aes(x=4.8, y=106, xend=5, yend=103), colour="darkgray")
mergedoffset <- mergedoffset + annotate("text", x=3.5, y=109, label="9.4%p (p=0.009)", size=3)
mergedoffset <- mergedoffset + scale_fill_manual(values=c("#CC0000","#FF0000","#CC33CC","#9900FF", "#0000FF"))
mergedoffset <- mergedoffset + theme(legend.position="none")
print(mergedoffset)
dev.off()



##### Table 2: ATEs for Linguistic Differences

# Hypothesis 1: Bystanders help immigrants speaking a foreign language less than natives
merged$hijab <- NA
merged$hijab[merged$languageoff==1 | merged$languageoff==2] <- 1
merged$hijab[merged$languageoff==3 | merged$languageoff==4 | merged$languageoff==5] <-0

merged$foreign <- NA
merged$foreign[merged$languageoff==1 | merged$languageoff==3] <- 1
merged$foreign[merged$languageoff==2 | merged$languageoff==4 | merged$languageoff==5] <- 0


ifvn   <- felm(help ~ hijab*foreign | 0 | 0 | 0, data=merged)
ifvnFE <- felm(help ~ hijab*foreign | exp + bystander | 0 | 0, data=merged)
ifvnFEt <- felm(help ~ hijab*foreign | exp + bystander + team | 0 | 0, data=merged)
summary(ifvn)
summary(ifvnFE)
summary(ifvnFEt)

# Hypothesis 2: Bystanders are less likely to help immigrants speaking a foreign language than immigrants speaking German

merged$ifvig <- NA
merged$ifvig[merged$languageoff==1 | merged$languageoff==3] <- 1
merged$ifvig[merged$languageoff==2 | merged$languageoff==4] <- 0

ifvig   <-  felm(help ~ ifvig | 0 | 0 | 0, data=merged)
ifvigFE <-  felm(help ~ ifvig | exp + bystander | 0 | 0, data=merged)
ifvigFEt <- felm(help ~ ifvig | exp + bystander + team | 0 | 0, data=merged)
summary(ifvig)
summary(ifvigFE)

stargazer(ifvn, ifvnFE, ifvnFEt, ifvig, ifvigFE, ifvigFEt,
          #dep.var.caption = c("Immigrant Foreign Language vs Native", "Immigrant Foreign Language vs Immigrant German "),
          dep.var.labels  = c("Any help?"),
          type="text", 
          column.separate = c(3,3),
          omit.stat = c("f"),
          keep=c("hijab", "foreign", "hijab:foreign", "ifvig",  "Constant"),
          keep.stat=c("n", "rsq"),
          add.lines = list(c("Sample", "Full", "Full", "Full", "Immigrant", "Immigrant", "Immigrant"),
                           c("Experiment FE", "No", "Yes", "Yes", "No", "Yes", "Yes"),
                           c("Bystander FE",  "No", "Yes", "Yes", "No", "Yes", "Yes"),
                           c( "Team FE", "No", "No", "Yes", "No", "No", "Yes")))
